In [ ]:
# importing libraries
 
import numpy as np
import matplotlib.pyplot as plt
from tqdm.auto import trange, tqdm
from matplotlib import cm
from matplotlib.colors import Normalize, ListedColormap, LinearSegmentedColormap
from matplotlib.cm import ScalarMappable
import random
from scipy.optimize import curve_fit
import time
import json

plt.style.use("default")
#plt.style.use("dark_background")
In [ ]:
%%html
<style>
.cell-output-ipywidget-background {
    background-color: transparent !important;
}
:root {
    --jp-widgets-color: var(--vscode-editor-foreground);
    --jp-widgets-font-size: var(--vscode-editor-font-size);
}  
</style>
In [ ]:
def infer_style():
    for style_name, style_params in plt.style.library.items():
        if all(key in plt.rcParams and plt.rcParams[key] == val for key, val in style_params.items()):
            return style_name
    return 'Default'
In [ ]:
# huge python class with everything inside
#
# m = (m_r, m_g, m_b) 3 colors
# the grid is from (0,0) to (n,n) make NxN square, each is a 2 dim list, each vertex has an horizontal and vertical edges like so:
#               0
#             ----(x,y)
#                   |    1
#                   |
#
# the 4 directions are encoded as:
#                     0
#                  -------
#             3    |     |   1
#                  |_____|
#                     2
#

class stateSpace:
    def __init__(self, num_colors, grid_size, beta, init = 0, bc = 0, algo = 'metropolis'):  # 'glauber'      
    
        self.grid_size = grid_size
        self.V = grid_size**2
        self.num_colors = num_colors
        self.beta = beta
        self.accepted = 0
        self.rejected = 0
        self.algo = algo
        
        self.data = {}  # dictionary to store data
        
        self.shape = (num_colors, grid_size+2, grid_size+2, 2)  #the +2 is for free bc  we work only on vertices from  (1,grid_size+1)

        self.grid = 2*np.random.randint(0, 10, self.shape, dtype=int) if bc == 'random' else bc*np.ones(self.shape, dtype=int)

        if init == 'random':
            self.random_init()
            
        else:
            self.uniform_init(init)

    def random_init(self):
        self.grid[:, 1:-1, 1:-1, :] = 2*np.random.randint(0, 10, size = (self.num_colors, self.grid_size, self.grid_size, 2), dtype=int)
        # bottom border
        self.grid[:, 1:-1, 0, 0] = 2*np.random.randint(0, 10, size = (self.num_colors, self.grid_size))
        # left border
        self.grid[:, 0, 1:-1, 1] = 2*np.random.randint(0, 10, size = (self.num_colors, self.grid_size))

    def uniform_init(self, k):
        self.grid[:, 1:-1, 1:-1, :] = k*np.ones((self.num_colors, self.grid_size, self.grid_size, 2), dtype=int)
        
    def step(self, num_steps = 1, progress_bar = True, sample_rate = 10_000, observables = []):   # observables is a list of functions like [m.avg_links, m.avg_local_time]
        
        # init the dictionary and add some info
        self.data['steps'] = num_steps
        self.data['beta'] = self.beta
        self.data['grid_size'] = self.grid_size
        
        if observables != []: 
            for ob in observables:
                self.data[ob.__name__] = [] # add keys to the dictionary for each observable
                
        for i in trange(num_steps, disable = not progress_bar):  
            # store data every sample_rate steps
            if i % sample_rate == 0:   
                for ob in observables:
                    self.data[ob.__name__].append(ob().tolist())

            # choose a random color 
            c = np.random.randint(0, self.num_colors, dtype=int)
         
            # choose a random square
            s = np.random.randint(1, self.grid_size+1, size = 2, dtype=int)

            # get num_link on each side of square s
            S = np.zeros(4, dtype=int)
            S[0] = self.grid[c, s[0], s[1], 0]
            S[2] = self.grid[c, s[0], s[1]-1, 0]
            S[1] = self.grid[c, s[0], s[1], 1]
            S[3] = self.grid[c, s[0]-1, s[1], 1]

            # get list of all possible transformation
            transformations = self.get_possible_transformations(S)                    ############################ MINIMAL OR FULL ST ####################################
            #transformations = self.minimal_transformations(S)
            
            # pick uniformly a random transformation
            M = len(transformations)   # num of possible transformation of current state, compute only once! we also need it to compute tha ratio M/M_prime in acceptance_prob
            index = np.random.randint(0, M)  
            X = transformations[index]
            #print('transformation: {}'.format(X))
            
            if self.acceptance_prob(S, M, s, X) >= random.random():
                self.accepted += 1
                self.square_transformation(c, s, X)
            else:
                self.rejected += 1
    
    def minimal_transformations(self, S):
        # list of just minimal transformations for irreducibility

        transformations = [ 
            (1, 1, 1, 1),                      # uniform +1
            (-2, 0, 0, 0),                    # single top
            (0, -2, 0, 0),                    # single right
            (0, 0, -2, 0),                    # single bottom
            (0, 0, 0, -2)                    # single left
            ]    
        
        if S[0] < 2:                                                # top is 0 or 1, remove single top -2
            transformations.remove( (-2, 0, 0, 0) )
            
        if S[1] < 2:
            transformations.remove((0, -2, 0, 0))
           
        if S[2] < 2:
            transformations.remove((0, 0, -2, 0))
        
        if S[3] < 2:
            transformations.remove((0, 0, 0, -2))

        return transformations 
    
    def get_possible_transformations(self, S):
        # list of all possible transformations

        transformations = [ 
            (-1,-1,-1,-1),                    # uniform                     
            (-2, 0, 0, 0),                    # single top
            (0, -2, 0, 0),                    # single right
            (0, 0, -2, 0),                    # single bottom
            (0, 0, 0, -2),                    # single left
            (1,-1, 1,-1), (-1, 1,-1, 1),      # swap opposite
            (1, 1,-1,-1), (-1,-1, 1, 1),      # swap adjacent
            (-1, 1,1,-1), ( 1,-1,-1, 1)#,   
            #(1,-1,-1,-1), (-1, 1, 1, 1),      # triple up
            #(-1,1,-1,-1), ( 1,-1, 1, 1),      # triple right
            #(-1,-1,1,-1), ( 1, 1,-1, 1),      # triple bottom
           # (-1,-1,-1,1), ( 1, 1, 1,-1)       # triple left
            ]    
        
        if S[0] < 2:                                                # top is 0 or 1, remove single top -2
            transformations.remove( (-2, 0, 0, 0) )
            if S[0] == 0:                                           # top is 0, remove uniform -1, swap 
                transformations.remove((-1, -1, -1, -1))
                transformations.remove((-1, 1, -1, 1))
                transformations.remove((-1, -1, 1, 1))
                transformations.remove((-1, 1, 1, -1))
                #transformations.remove((-1, 1, 1, 1))
                #transformations.remove((-1, 1, -1, -1))
                #transformations.remove((-1, -1, 1, -1))
                #transformations.remove((-1, -1, -1, 1))
        if S[1] < 2:
            transformations.remove((0, -2, 0, 0))
            if S[1] == 0:
                if (-1, -1,-1,-1) in transformations: transformations.remove((-1, -1,-1,-1))
                if (1, -1, 1, -1) in transformations: transformations.remove(( 1, -1, 1,-1))
                if (-1, -1, 1, 1) in transformations: transformations.remove((-1, -1, 1, 1))
                if (1, -1, -1, 1) in transformations: transformations.remove(( 1, -1,-1, 1))
                #if (1, -1, -1, -1) in transformations: transformations.remove((1, -1, -1, -1))
                #if (1, -1, 1, 1) in transformations: transformations.remove((1, -1, 1, 1))
                #if (-1, -1, 1, -1) in transformations: transformations.remove((-1, -1, 1, -1))
                #if (-1, -1, -1, 1) in transformations: transformations.remove((-1, -1, -1, 1))

        if S[2] < 2:
            transformations.remove((0, 0, -2, 0))
            if S[2] == 0:
                if (-1, -1,-1,-1) in transformations: transformations.remove((-1,-1, -1,-1))
                if (-1, 1, -1, 1) in transformations: transformations.remove((-1, 1, -1, 1))
                if (1, 1, -1, -1) in transformations: transformations.remove(( 1, 1, -1,-1))
                if (1, -1, -1, 1) in transformations: transformations.remove(( 1,-1, -1, 1))
                #if (1, -1, -1, -1) in transformations: transformations.remove((1, -1, -1, -1))
                #if (-1, 1, -1, -1) in transformations: transformations.remove((-1, 1, -1, -1))
                #if (1, 1, -1, 1) in transformations: transformations.remove((1, 1, -1, 1))
                #if (-1, -1, -1, 1) in transformations: transformations.remove((-1, -1, -1, 1))
        
        if S[3] < 2:
            transformations.remove((0, 0, 0, -2))
            if S[3] == 0:
                if (-1, -1,-1,-1) in transformations: transformations.remove((-1,-1,-1, -1))
                if (1, -1, 1, -1) in transformations: transformations.remove(( 1,-1, 1, -1))
                if (1, 1, -1, -1) in transformations: transformations.remove(( 1, 1,-1, -1))
                if (-1, 1, 1, -1) in transformations: transformations.remove((-1, 1, 1, -1))
                #if (1, -1, -1, -1) in transformations: transformations.remove((1, -1, -1, -1))
                #if (-1, 1, -1, -1) in transformations: transformations.remove((-1, 1, -1, -1))
                #if (-1, -1, 1, -1) in transformations: transformations.remove((-1, -1, 1, -1))
                #if (1, 1, 1, -1) in transformations: transformations.remove((1, 1, 1, -1))


        return transformations + [(1, 1, 1, 1), (2, 0, 0, 0), (0, 2, 0, 0), (0, 0, 2, 0), (0, 0, 0, 2)]   #add back always applicable transformations

    def square_transformation(self, c, s, X):    # X = (a_1, a_2, a_3, a_4) like in the thesis
        #top
        self.grid[c, s[0], s[1], 0] += X[0]
        #right
        self.grid[c, s[0], s[1], 1] += X[1]
        #bottom
        self.grid[c, s[0], s[1]-1, 0] += X[2]
        #left
        self.grid[c, s[0]-1, s[1], 1] += X[3]
    
    def acceptance_prob(self, S, M, s, X): # S = (l1,l2,l3,l4) links on the square, M, s = (x,y) vertex in the grid, X = (a_1,a_2,a_3,a_4) square transformation
        # possible transformation ratio
        S_prime = np.copy(S) 
        S_prime += np.array(X) #apply the transformation in the square

        #get M_prime, the number of possibile transformation of the new state
        M_prime = len(self.get_possible_transformations(S_prime))                                 ############################ MINIMAL OR FULL ST ####################################
        #M_prime = len(self.minimal_transformations(S_prime))

        # prob ratio
        match X:
            case (1, 1, 1, 1):
                A = self.beta**4 / (16 * (S[0] + 1)*(S[1] + 1)*(S[2] + 1)*(S[3] + 1)*(self.num_colors/2 + self.get_local_time(s[0], s[1]))*(self.num_colors/2 + self.get_local_time(s[0], s[1]-1))*(self.num_colors/2 + self.get_local_time(s[0]-1 , s[1]))*(self.num_colors/2 + self.get_local_time(s[0]-1 , s[1]-1)))
            case (-1,-1,-1,-1):
                A = (16 / (self.beta**4)) * S[0]*S[1]*S[2]*S[3] *(self.num_colors/2 + self.get_local_time(s[0], s[1]) - 1)*(self.num_colors/2 + self.get_local_time(s[0], s[1]-1 ) - 1)*(self.num_colors/2 + self.get_local_time(s[0]-1, s[1]) - 1)*(self.num_colors/2 + self.get_local_time(s[0]-1, s[1]-1) - 1)
            case (2, 0, 0, 0):
                A = self.beta**2/ (4*(S[0]+2)*(S[0]+1)*(self.num_colors/2 + self.get_local_time(s[0], s[1]))*(self.num_colors/2 + self.get_local_time( s[0]-1, s[1]))  )
            case (-2, 0, 0, 0):
                A = 4 * S[0]*(S[0]-1) / self.beta**2 * (self.num_colors/2 + self.get_local_time(s[0], s[1]) - 1)*(self.num_colors/2 + self.get_local_time(s[0]-1, s[1]) - 1)
            case (0, 2, 0, 0):
                A = self.beta**2/ (4*(S[1]+2)*(S[1]+1)*(self.num_colors/2 + self.get_local_time(s[0], s[1]))*(self.num_colors/2 + self.get_local_time(s[0], s[1]-1))  )
            case (0,-2, 0, 0): 
                A = 4 * S[1]*(S[1]-1) / self.beta**2 * (self.num_colors/2 + self.get_local_time(s[0], s[1]) - 1)*(self.num_colors/2 + self.get_local_time(s[0], s[1]-1) - 1)
            case (0, 0, 2, 0):
                A = self.beta**2/ (4*(S[2]+2)*(S[2]+1)*(self.num_colors/2 + self.get_local_time( s[0]-1, s[1]-1))*(self.num_colors/2 + self.get_local_time(s[0], s[1]-1))  )
            case (0, 0,-2, 0):
                A = 4 * S[2]*(S[2]-1) / self.beta**2 * (self.num_colors/2 + self.get_local_time( s[0]-1, s[1]-1 ) - 1)*(self.num_colors/2 + self.get_local_time(s[0], s[1]-1 ) - 1)
            case (0, 0, 0, 2):
                A = self.beta**2/ (4*(S[3]+2)*(S[3]+1)*(self.num_colors/2 + self.get_local_time( s[0]-1, s[1]))*(self.num_colors/2 + self.get_local_time( s[0]-1, s[1]-1))  )
            case (0, 0, 0,-2):
                A = 4 * S[3]*(S[3]-1) / self.beta**2 * (self.num_colors/2 + self.get_local_time( s[0]-1, s[1]) - 1)*(self.num_colors/2 + self.get_local_time( s[0]-1, s[1]-1) - 1)
            
            case (-1, 1, -1, 1):
                A = S[0]*S[2]/( (S[1]+1)*(S[3]+1) )
            case (1, -1, 1, -1):
                A = S[1]*S[3]/( (S[0]+1)*(S[2]+1) )
                
            case (-1, -1, 1, 1):
                A = S[0]*S[1] / ( (S[2]+1)*(S[3]+1) ) * (self.num_colors/2 + self.get_local_time(s[0], s[1]) - 1 ) / (self.num_colors/2 + self.get_local_time(s[0]-1, s[1]-1))
            case (1, 1, -1, -1):
                A = S[2]*S[3] / ( (S[0]+1)*(S[1]+1) ) * (self.num_colors/2 + self.get_local_time(s[0]-1, s[1]-1) -1 ) / (self.num_colors/2 + self.get_local_time(s[0], s[1])) 
            case (-1, 1, 1, -1):
                A = S[0]*S[3] / ( (S[1]+1)*(S[2]+1) ) * (self.num_colors/2 + self.get_local_time(s[0]-1, s[1]) -1 ) / (self.num_colors/2 + self.get_local_time(s[0], s[1]-1)) 
            case (1, -1, -1, 1):
                A = S[2]*S[1] / ( (S[3]+1)*(S[0]+1) ) * (self.num_colors/2 + self.get_local_time(s[0], s[1]-1) -1 ) / (self.num_colors/2 + self.get_local_time(s[0]-1, s[1]) ) 

            # missing triple transformations

        #print('acceptance prob = {}'.format(min(1, M/M_prime * A)))
        return min(1, M/M_prime * A) if self.algo == 'metropolis' else 1/(1 + M_prime/(M*A))   # Metropolis  Glauber       #### May impact performanca a bit, better to edit the code!
    
    def get_local_time(self, x, y):   # we know already the number of links in square s! we are wasting a bit of compute power 
        local_time = 0
        for c in range(self.num_colors):
            local_time += self.grid[c, x, y, 0] + self.grid[c, x, y, 1] + self.grid[c, x, y + 1, 1] + self.grid[c, x + 1, y, 0]
        return local_time // 2

    def max_links(self):
        max_links = np.zeros(self.num_colors)
        for c in range(0,self.num_colors):
            for x in range(1,self.grid_size+1):
                for y in range(1,self.grid_size+1):
                    if self.grid[c, x, y, 0] >= max_links[c]:
                        max_links[c] = self.grid[c, x, y, 0]
                    if self.grid[c, x, y, 1] >= max_links[c]:
                        max_links[c] = self.grid[c, x, y, 1]
        return max_links
    
    def avg_links(self):
        return np.mean(self.grid[:, 1:-1, 1:-1, :])

    def avg_local_time(self):
        total_local_time = 0
        for x in range(1,self.grid_size+1):
            for y in range(1,self.grid_size+1):
                total_local_time += self.get_local_time(x,y)
        return total_local_time / self.V
    
    def loop_builder(self):  # given a link configuration, randomly builds loops (we are choosing a link pairing uniformly) returns list of loops (as a sequence of vertices) for each color
        loops = []
        lenghts = []
        for c in range(self.num_colors):
            color_loops = []
            #copy the grid 
            G = np.copy(self.grid[c])
            nz = G.nonzero()
            non_zero = len(nz[0])
            if non_zero == 0:
                    return [], [0]
            
            while non_zero > 0:
                #pick first non-zero and unvisited edge
                nz = G.nonzero()
                non_zero = len(nz[0])
                if non_zero == 0:
                    break
                # choose a random vertex with non-zero incident links 
                rand_index = np.random.randint(0, non_zero)
                x, y  = nz[0][rand_index], nz[1][rand_index]
                starting_vertex = (x,y)
                current_loop = []
                
                # first step outside loops
                dir = []
                
                top = G[x,y+1,1]
                right = G[x+1,y,0]
                bottom = G[x,y,1]
                left = G[x,y,0]
                
                if top > 0:
                    dir.extend([0]*top)
                elif right > 0:
                    dir.extend([1]*right)
                elif bottom > 0:
                    dir.extend([2]*bottom)
                elif left > 0:
                    dir.extend([3]*left)
                
                # pick a random dir with prob prop to num_links  
                rand_dir = np.random.choice(dir)
                
                match rand_dir:
                    case 0:
                        # remove one link
                        G[x,y+1,1] -= 1
                        #move there
                        y += 1
                    case 1:
                        G[x+1,y,0] -= 1
                        x += 1
                    case 2:
                        G[x,y,1] -= 1
                        y -= 1
                    case 3:
                        G[x,y,0] -= 1
                        x -= 1
                length = 0
                while True:
                    current_loop.append((x,y))
                    length += 1
                    # look if we can trav in each of the 4 directions: top = 0, right = 1, down = 2 and left = 3 with prob eq. to num_links/Z
                    dir = []
                    
                    top = G[x,y+1,1]
                    right = G[x+1,y,0]
                    bottom = G[x,y,1]
                    left = G[x,y,0]
                    
                    if top > 0:
                        dir.extend([0]*top)
                    elif right > 0:
                        dir.extend([1]*right)
                    elif bottom > 0:
                        dir.extend([2]*bottom)
                    elif left > 0:
                        dir.extend([3]*left)

                    if (x,y) == starting_vertex:
                        if random.random() <= 1/(len(dir)+1):
                            lenghts.append(length)
                            break
                    
                    # pick a random dir with prob prop to num_links  
                    rand_dir = np.random.choice(dir)
                    
                    match rand_dir:
                        case 0:
                            # remove one link
                            G[x,y+1,1] -= 1
                            #move there
                            y += 1
                        case 1:
                            G[x+1,y,0] -= 1
                            x += 1
                        case 2:
                            G[x,y,1] -= 1
                            y -= 1
                        case 3:
                            G[x,y,0] -= 1
                            x -= 1
                        
                color_loops.append(current_loop)
            loops.append(color_loops)
        return loops, lenghts
    
    def avg_loop_length(self):
        _, lengths = self.loop_builder()
        return np.mean(lengths)
                
    def check_state(self): # checks if the current state m is legal (every vertex has even degree)
        for c in range(self.num_colors):
            for x in range(1,self.grid_size+1):
                for y in range(1,self.grid_size+1):
                    if (self.grid[c, x, y, 0] + self.grid[c, x, y, 1] + self.grid[c, x, y + 1, 1] + self.grid[c, x + 1, y, 0]) % 2 != 0:
                        print('###  illegal state!  ###')
                        return False 
        return True 
    def plot_one_color(self, c, cmap, ax, alpha = 1.0): # plots the grid of just one color
        for x in range(0, self.grid_size+2):
                for y in range(0, self.grid_size+2):
                    # horizontal
                    if self.grid[c,x,y,0] != 0:
                        edge_color = cmap(self.grid[c,x,y,0])
                        ax.plot([x-1, x], [y, y], color=edge_color,linewidth=1.5, alpha = alpha)
                    # vertical
                    if self.grid[c,x,y,1] != 0:   
                        edge_color = cmap(self.grid[c,x,y,1])
                        ax.plot([x, x], [y, y-1], color=edge_color, linewidth=1.5, alpha = alpha)
        
    def plot_loop(self, c, loop, color = 'yellow', alpha = 0.25):  # plots the longest loop over the grid
        fig, ax = plt.subplots(figsize=(12,12))
        num_segments = int(self.max_links()[c]+1)            #color dependet!
        cmap = create_cmap(c, num_segments)
        self.plot_one_color(c, cmap, ax)
        for i in range(len(loop)-1):
            ax.plot( [loop[i][0], loop[i+1][0]], [loop[i][1], loop[i+1][1]], linewidth=1.5, color = color, alpha = alpha)
        #draw last link
        ax.plot( [loop[0][0], loop[-1][0]], [loop[0][1], loop[-1][1]], linewidth=1.5, color = color, alpha = alpha)
        ax.set_title('length = {}'.format(len(loop)))
        
    def plot_grid(self, figsize=(32, 8), colorbar = True, file_name = None):     # plots the grid of every color 
        fig, axes = plt.subplots(1,self.num_colors,figsize = figsize, gridspec_kw={'hspace': 0.05, 'wspace': 0.05}) #, facecolor='black')
        # Adjust the space between subplots
        for c in range(self.num_colors):
            # Define a colormap
            num_segments = int(self.max_links()[c]+1)            #color dependet!
            cmap = create_cmap(c, num_segments)
            
            # Create a ScalarMappable for colorbar
            norm = Normalize(vmin=0, vmax=num_segments)
            sm = ScalarMappable(cmap=cmap, norm=norm)
            sm.set_array([])  # Empty array, as we'll not use actual data

            self.plot_one_color(c, cmap, axes[c])
            #axes[c].set_title('avg links = {}'.format(self.avg_links()[c]))
            axes[c].set_xlim(-(1+0.05*self.grid_size), 2+self.grid_size*1.05)
            axes[c].set_ylim(-(1+0.05*self.grid_size), 2+self.grid_size*1.05)
            #axes[c].axis('square')
            axes[c].axis('off')                                                ########### axis
            
            # Add colorbar
            if colorbar:
                # Add colorbar
                cbar = plt.colorbar(sm, ax=axes[c])
                cbar.set_ticks(  0.5 + np.arange(0, num_segments,1))
                cbar.set_ticklabels(list(range(0, num_segments)))
                #cbar.set_label('Color Mapping')

        fig.suptitle(r'grid size = {}     $\beta$ = {}        steps = {:g}'.format(self.grid_size, self.beta, self.accepted + self.rejected))
        #save it
        if file_name != None:
            plt.savefig(file_name)
        plt.show()

    def plot_overlap(self, figsize = (12,12), normalized = False, file_name = None):  # plots every color overlapped
        # Create a figure and axes
        fig, ax = plt.subplots(figsize = figsize)

        for c in range(self.num_colors):
            # Define a colormap
            num_segments = int(self.max_links()[c]+1) if not normalized else 2
            cmap = create_cmap(c, num_segments)
            
            # Create a ScalarMappable for colorbar
            norm = Normalize(vmin=0, vmax=num_segments)
            sm = ScalarMappable(cmap=cmap, norm=norm)
            sm.set_array([])  # Empty array, as we'll not use actual data

            self.plot_one_color(c, cmap, ax, 0.6)
            ax.set_title(r'grid size = {}     $\beta$ = {}        steps = {:g}'.format(self.grid_size, self.beta, self.accepted + self.rejected))
            ax.set_xlim(-(1+0.05*self.grid_size), 2+self.grid_size*1.05)
            ax.set_ylim(-(1+0.05*self.grid_size), 2+self.grid_size*1.05)
            
            #ax.axis('square')
            
            ax.set_xticks(np.arange(1,self.grid_size+1), minor = True)
            ax.set_yticks(np.arange(1,self.grid_size+1), minor = True)
            ax.grid(which='both')
            #ax.axis('off')
         
            
    
        #save it
        if file_name != None:
            plt.savefig(file_name)
        plt.show()
    
    def summary(self):   # prints some stats
        print('average number of links: {}'.format(self.avg_links()))
        print('max number of links: {}'.format(max(self.max_links()) ))
        print('avg local time: {}'.format(self.avg_local_time()))
        loops = self.loop_builder()[0] #stats on color 0
        lenghts = [len(l) for l in loops]
        print('avg loop length: {}'.format(np.mean(lenghts)))
        print('max loop length: {}'.format(max(lenghts)))
        steps = self.accepted + self.rejected
        if steps == 0:
            print('steps = {:g}   acceptance ratio = {:.6f}'.format(steps, 0))
        else:
            print('steps = {:g}   acceptance ratio = {:.6f}'.format(steps, self.accepted / (steps)))
    
def create_cmap(color, n_bins):  # create a color map to show how many links are in a edge
    # get current style
    current_style = infer_style()
    start = (0,0,0) if current_style == 'dark_background' else (1,1,1)
    if color == 0:
        colors = [start, (1, 0, 0)]  # RGB values  
    if color == 1:
        colors = [start, (0, 1, 0)]  # RGB values
    if color == 2:
        colors = [start, (0, 0, 1)]  # RGB values

    cmap = LinearSegmentedColormap.from_list('my_cmap', colors, N=n_bins)
    return cmap
In [ ]:
file_name = 'data.json'

with open(file_name, 'w') as file:
    json.dump(m.data, file)
In [ ]:
# initialize the grid
m = stateSpace(num_colors = 3, grid_size = 32, beta = 5 , init = 0, bc = 0, algo='metropolis')  

# run it
m.step(100_000, sample_rate= 10_000, observables = [m.avg_links, m.avg_local_time, m.avg_loop_length])  # about 3 seconds for 100k steps on my PC

# plot it
m.plot_grid() #, file_name = '64_0_9_U2.pdf')

#sample loops
loops = m.loop_builder()

m.summary()
  0%|          | 0/100000 [00:00<?, ?it/s]
No description has been provided for this image
average number of links: 0.5457356770833334
max number of links: 5.0
avg local time: 3.2294921875
avg loop length: 179.66666666666666
max loop length: 192
steps = 100000   acceptance ratio = 0.237620
In [ ]:
m.data
Out[ ]:
{'steps': 100000,
 'beta': 5,
 'grid_size': 32,
 'avg_links': [0.0,
  0.4482421875,
  0.4957682291666667,
  0.5152994791666666,
  0.5322265625,
  0.53857421875,
  0.5465494791666666,
  0.5343424479166666,
  0.5328776041666666,
  0.55126953125],
 'avg_local_time': [0.0,
  2.646484375,
  2.9296875,
  3.0361328125,
  3.1416015625,
  3.177734375,
  3.224609375,
  3.1552734375,
  3.146484375,
  3.2578125],
 'avg_loop_length': [0.0,
  3.7244094488188977,
  4.271370420624152,
  5.153238546603475,
  5.722033898305085,
  5.634266886326194,
  5.83445945945946,
  6.416666666666667,
  6.337711069418386,
  6.087108013937282]}

Equilibrium check¶

We sample every $10k$ steps the average number of links, average local time and average loop length. We plot these values normalized as a functions of time to check if the chain has reached equilibrium.

In [ ]:
%%time

STEPS = 1_000_000
SAMPLE_RATE = 10_000

# define an array of betas
betas = [0.5, 1, 4, 8, 32, 64]

# plot relevant observables normalized to last value
fig, axes = plt.subplots(3, 2, figsize = (16,14))

for i in range(3):
    for j in range(2):
        # initialize the grid and run the simulation
        m = stateSpace(num_colors = 3, grid_size = 32, beta = betas[i+3*j], init = 0, bc = 0, algo = 'metropolis')  
        observables = [m.avg_links, m.avg_local_time, m.avg_loop_length]
        m.step(num_steps = STEPS, progress_bar = False, sample_rate = SAMPLE_RATE, observables = observables)  
        
        x = np.arange(0, STEPS, SAMPLE_RATE)  
        
        # normalize with avg of last 10 values 
        for ob in observables:
            avg = np.mean(m.data[ob.__name__][-10:-1])
            axes[i,j].plot(x, np.array(m.data[ob.__name__]) / avg , label = ob.__name__)
            
        axes[i,j].legend()
        axes[i,j].set_xlabel('steps')
        axes[i,j].set_title(r'$\beta$ = {}'.format(m.beta))

        axes[i,j].grid(linewidth=0.5, linestyle = '--')
        axes[i,j].set_xticks(np.arange(0, STEPS, STEPS // 10))
        axes[i,j].ticklabel_format(axis = 'x', style = 'sci', scilimits=(0,0))
        
fig.suptitle(r'Normalized observables       grid size = {}'.format(m.grid_size))
plt.savefig('equilibrium_check.pdf')
plt.show()
No description has been provided for this image
CPU times: total: 2min 53s
Wall time: 3min 9s
In [ ]:
print(m.data)
{'steps': 100000, 'beta': 5, 'grid_size': 32, 'avg_links': [0.0, 0.4482421875, 0.4957682291666667, 0.5152994791666666, 0.5322265625, 0.53857421875, 0.5465494791666666, 0.5343424479166666, 0.5328776041666666, 0.55126953125], 'avg_local_time': [0.0, 2.646484375, 2.9296875, 3.0361328125, 3.1416015625, 3.177734375, 3.224609375, 3.1552734375, 3.146484375, 3.2578125], 'avg_loop_length': [0.0, 3.7244094488188977, 4.271370420624152, 5.153238546603475, 5.722033898305085, 5.634266886326194, 5.83445945945946, 6.416666666666667, 6.337711069418386, 6.087108013937282]}

Long sim¶

We now focus on small and big values of $\beta$, namely 2 and 64, but with $10M$ steps.

In [ ]:
# low beta

# initialize the grid
m = stateSpace(num_colors = 3, grid_size = 64, beta = 2 , init = 0, bc = 0, algo='metropolis')  

# run it
m.step(10_000_000, sample_rate = 100_000)  # about 3 seconds for 100k steps on my PC

# plot relevant observables normalized to last value
fig, ax = plt.subplots(figsize = (16,8))

# steps
x = np.arange(0, 10_000_000, 100_000)  # we sampled every 10k steps

# normalize with avg of last 10 values 
n1 = np.mean(m.data1[-10:-1])
n2 = np.mean(m.data2[-10:-1])
n3 = np.mean(m.data3[-10:-1])


ax.plot(x, np.array(m.data1) / n1 , label = 'avg links')
ax.plot(x, np.array(m.data2) / n2, label = 'avg local time')
ax.plot(x, np.array(m.data3) / n3, label = 'avg loop length')
ax.legend()
ax.set_xlabel('steps')
ax.set_title(r'normalized observables              grid size = {}   $\beta$ = {}'.format(m.grid_size, m.beta))

ax.grid(linewidth=0.5, linestyle = '--')
ax.set_xticks(np.arange(0,10_100_000, 1_000_000))
ax.ticklabel_format(axis = 'x', style = 'sci', scilimits=(0,0))

plt.show()
  0%|          | 0/10000000 [00:00<?, ?it/s]
No description has been provided for this image
In [ ]:
# high beta

# initialize the grid
m = stateSpace(num_colors = 3, grid_size = 64, beta = 64 , init = 0, bc = 0, algo = 'metropolis')  

# run it
m.step(10_000_000, sample_rate = 100_000)  # about 3 seconds for 100k steps on my PC

# plot relevant observables normalized to last value
fig, ax = plt.subplots(figsize = (16,8))

# steps
x = np.arange(0, 10_000_000, 100_000)  # we sampled every 10k steps

# normalize with avg of last 10 values 
n1 = np.mean(m.data1[-10:-1])
n2 = np.mean(m.data2[-10:-1])
n3 = np.mean(m.data3[-10:-1])


ax.plot(x, np.array(m.data1) / n1 , label = 'avg links')
ax.plot(x, np.array(m.data2) / n2, label = 'avg local time')
ax.plot(x, np.array(m.data3) / n3, label = 'avg loop length')
ax.legend()
ax.set_xlabel('steps')
ax.set_title(r'normalized observables              grid size = {}   $\beta$ = {}'.format(m.grid_size, m.beta))

ax.grid(linewidth=0.5, linestyle = '--')
ax.set_xticks(np.arange(0,10_100_000, 1_000_000))
ax.ticklabel_format(axis = 'x', style = 'sci', scilimits=(0,0))

plt.show()
  0%|          | 0/10000000 [00:00<?, ?it/s]
No description has been provided for this image

Hand-crafted starting configurations¶

We now consider two hand-crafted initial conditions: one in which there is a loop that visits all edges, which we call snake, and one with uniform links but with an hole in the middle, which we call donut.

In [ ]:
STEPS = 1_000_000
SAMPLE_RATE = 10_000

m = stateSpace(num_colors = 3, grid_size = 32, beta = 8 , init = 0, bc = 0, algo = 'metropolis')

#building the donut state

for c in range(m.num_colors):
    for x in range(1, m.grid_size // 4+1):
        for y in range(1, m.grid_size+1):
            m.grid[c, x, y, 0] = 1
            m.grid[c, x, y, 1] = 1
    for x in range(m.grid_size+1 - m.grid_size // 4, m.grid_size+1):
        for y in range(1, m.grid_size+1):
            m.grid[c, x, y, 0] = 1
            m.grid[c, x, y, 1] = 1
for c in range(m.num_colors):
    for y in range(1, m.grid_size // 4+1):
        for x in range(1, m.grid_size+1):
            m.grid[c, x, y, 0] = 1
            m.grid[c, x, y, 1] = 1
    for y in range(m.grid_size+1 - m.grid_size // 4, m.grid_size + 1):
        for x in range(1, m.grid_size+1):
            m.grid[c, x, y, 0] = 1
            m.grid[c, x, y, 1] = 1

# add missing links on the outside border
for c in range(m.num_colors):
    for y in range(1, m.grid_size+1):
        m.grid[c,0, y, 1] = 1
        m.grid[c, y, 0, 0] = 1 
        
# add missing links on the inside border
for c in range(m.num_colors):
    for x in range(m.grid_size // 4, 3 * m.grid_size // 4 + 1):
        m.grid[c,x,3 * m.grid_size // 4, 0] = 1
        m.grid[c, 3 * m.grid_size // 4, x, 1] = 1 
        
#make it a legal state
m.grid *= 2

m.plot_overlap()
        
# run it
m.step(STEPS, sample_rate = SAMPLE_RATE)  # about 3 seconds for 100k steps on my PC

# plot relevant observables normalized to last value
fig, ax = plt.subplots(figsize = (16,8))

# steps
x = np.arange(0, 1_000_000, 10_000)  # we sampled every 10k steps

# normalize with avg of last 10 values 
n1 = np.mean(m.data1[-10:-1])
n2 = np.mean(m.data2[-10:-1])
n3 = np.mean(m.data3[-10:-1])


ax.plot(x, np.array(m.data1) / n1 , label = 'avg links')
ax.plot(x, np.array(m.data2) / n2, label = 'avg local time')
ax.plot(x[1:], np.array(m.data3)[1:] / n3, label = 'avg loop length') # we remove the starting loop lenght
ax.legend()
ax.set_xlabel('steps')
ax.set_title(r'normalized observables              grid size = {}   $\beta$ = {}'.format(m.grid_size, m.beta))

ax.grid(linewidth=0.5, linestyle = '--')
ax.set_xticks(np.arange(0,STEPS, 10*SAMPLE_RATE))
ax.ticklabel_format(axis = 'x', style = 'sci', scilimits=(0,0))

plt.show()
No description has been provided for this image
  0%|          | 0/1000000 [00:00<?, ?it/s]
No description has been provided for this image
In [ ]:
m.plot_grid()
No description has been provided for this image
In [ ]:
m = stateSpace(num_colors = 3, grid_size = 32, beta = 8 , init = 0, bc = 0, algo = 'metropolis')

# building the snake state

# build concentric squares
offset = m.grid_size // 2+1
for c in range(m.num_colors):
    for k in range(1,offset):
        m.grid[c, k:-k, k-1, 0] = 1
        m.grid[c, k-1, k:-k, 1] = 1
        
        m.grid[c, (offset-k):(k-offset), offset + k -1, 0] = 1
        m.grid[c, offset + k -1, (offset-k):(k-offset), 1] = 1
        
# swap some links to join the squares
for c in range(m.num_colors):
    for k in range(1,offset-1):
        m.grid[c, k+1 , m.grid_size-k+1, 0] = 0
        m.grid[c, k+1 , m.grid_size-k, 0] = 0
        
        m.grid[c, k , m.grid_size-k+1, 1] = 1
        m.grid[c, k+1 , m.grid_size-k+1, 1] = 1
        
        

m.plot_overlap()
m.summary()

# run it
m.step(STEPS, sample_rate = SAMPLE_RATE)  # about 3 seconds for 100k steps on my PC

# plot relevant observables normalized to last value
fig, ax = plt.subplots(figsize = (16,8))

# steps
x = np.arange(0, 1_000_000, 10_000)  # we sampled every 10k steps

# normalize with avg of last 10 values 
n1 = np.mean(m.data1[-10:-1])
n2 = np.mean(m.data2[-10:-1])
n3 = np.mean(m.data3[-10:-1])


ax.plot(x, np.array(m.data1) / n1 , label = 'avg links')
ax.plot(x, np.array(m.data2) / n2, label = 'avg local time')
ax.plot(x[1:], np.array(m.data3)[1:] / n3, label = 'avg loop length')  # we remove the first avg loop lenght
ax.legend()
ax.set_xlabel('steps')
ax.set_title(r'normalized observables              grid size = {}   $\beta$ = {}'.format(m.grid_size, m.beta))

ax.grid(linewidth=0.5, linestyle = '--')
ax.set_xticks(np.arange(0,STEPS+ 10*SAMPLE_RATE, 10*SAMPLE_RATE))
ax.ticklabel_format(axis = 'x', style = 'sci', scilimits=(0,0))

plt.show()
No description has been provided for this image
average number of links: 0.5
max number of links: 1.0
avg local time: 2.9970703125
avg loop length: 1088.0
max loop length: 1088
steps = 0   acceptance ratio = 0.000000
  0%|          | 0/1000000 [00:00<?, ?it/s]
No description has been provided for this image
In [ ]:
m.plot_grid()
No description has been provided for this image

Loops¶

In [ ]:
# choose a color
color = 0

loops, lenghts = m.loop_builder()

plt.hist(lenghts, log = True, align='left', bins = range(min(lenghts), max(lenghts) + 2, 1))
plt.title('loop length distribution')
np.mean(lenghts)
Out[ ]:
6.476547842401501
No description has been provided for this image
In [ ]:
# find longest loop
max_length = max(lenghts)

longest_loop =  []

for l in loops[color]:
    if len(l) == max_length:
        longest_loop = l

#plot it 
m.plot_loop(color, longest_loop)
No description has been provided for this image

Sampling¶

In [ ]:
beta = 3
num_samples = 10
steps = 100_000
lenghts = []
mean_lengths = []
max_lengths = []
local_times = []

m = stateSpace(num_colors = 3, grid_size = 32, beta = beta , init = 0, bc = 0, algo='metropolis')  #(num_colors = 3, grid_size = 10, beta = 1, uniform_init = False)  #0.1875
m.step(100_000)

for i in trange(num_samples):
   
    m.step(steps, progress_bar=False)
    
    #loops
    loops = m.path_builder()

    # choose a color
    color = 0

    # compute lenghts
    lens = [len(l) for l in loops[color]]

    # find longest loop
    max_lengths.append(max(lens))
    lenghts.extend(lens)
    mean_lengths.append(np.mean(lens))
    
  0%|          | 0/100000 [00:00<?, ?it/s]
  0%|          | 0/10 [00:00<?, ?it/s]
In [ ]:
plt.hist(lenghts, log = True, density = True)
Out[ ]:
(array([0.06895581, 0.00211089, 0.00070363, 0.00105545, 0.        ,
        0.00035182, 0.        , 0.        , 0.        , 0.00035182]),
 array([  2. ,  15.6,  29.2,  42.8,  56.4,  70. ,  83.6,  97.2, 110.8,
        124.4, 138. ]),
 <BarContainer object of 10 artists>)
No description has been provided for this image

Animations¶

In [ ]:
def plot_one(grid, c, cmap, ax, alpha = 1.0):
        for x in range(0, len(grid[0])):
                for y in range(0, len(grid[0])):
                    # horizontal
                    if grid[c,x,y,0] != 0:
                        edge_color = cmap(grid[c,x,y,0])
                        ax.plot([x-1, x], [y, y], color=edge_color,linewidth=1.5, alpha = alpha)
                    # vertical
                    if grid[c,x,y,1] != 0:   
                        edge_color = cmap(grid[c,x,y,1])
                        ax.plot([x, x], [y, y-1], color=edge_color, linewidth=1.5, alpha = alpha)


def plot_over(i):
        figsize = (12,12) 
        normalized = False
        file_name = None
        grid_size = len(m.data1[0][0])
        # Create a figure and axes
        fig, ax = plt.subplots(figsize = figsize)

        for c in range(3):
            # Define a colormap
            num_segments = int(np.max(m.data1[i][c]) + 1) if not normalized else 2
            cmap = create_cmap(c, num_segments)
            
            # Create a ScalarMappable for colorbar
            norm = Normalize(vmin=0, vmax=num_segments)
            sm = ScalarMappable(cmap=cmap, norm=norm)
            sm.set_array([])  # Empty array, as we'll not use actual data

            plot_one(m.data1[i], c, cmap, ax, 0.6)
            #ax.set_title(r'grid size = {}     $\beta$ = {}        steps = {:g}'.format(self.grid_size, self.beta, self.accepted + self.rejected))
            ax.set_xlim(-(1+0.05*grid_size), 2+grid_size*1.05)
            ax.set_ylim(-(1+0.05*grid_size), 2+grid_size*1.05)
            #ax.axis('square')
            ax.set_title('step {}'.format( int(i*10_000)) )
            ax.axis('off')
    
        #save it
        if file_name != None:
            plt.savefig(file_name)
        plt.show()
'''
def anim_plotter(i):
    fig, ax = plt.subplots(figsize = (14,10)) 
    # Adjust the space between subplots

    # Define a colormap
    num_segments = np.max(m.data1[i][0])            #color dependet!
    cmap = create_cmap(0, num_segments)

    # Create a ScalarMappable for colorbar
    norm = Normalize(vmin=0, vmax=num_segments)
    sm = ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array([])  # Empty array, as we'll not use actual data

    grid_size = len(m.data1[0][0])

    plot_one(m.data1[i], 0, cmap, ax)
    #axes[c].set_title('avg links = {}'.format(self.avg_links()[c]))
    ax.set_xlim(-(1+0.05*grid_size), 2+grid_size*1.05)
    ax.set_ylim(-(1+0.05*grid_size), 2+grid_size*1.05)
    
    #axes[c].axis('square')
    #plt.show()
'''
Out[ ]:
"\ndef anim_plotter(i):\n    fig, ax = plt.subplots(figsize = (14,10)) \n    # Adjust the space between subplots\n\n    # Define a colormap\n    num_segments = np.max(m.data1[i][0])            #color dependet!\n    cmap = create_cmap(0, num_segments)\n\n    # Create a ScalarMappable for colorbar\n    norm = Normalize(vmin=0, vmax=num_segments)\n    sm = ScalarMappable(cmap=cmap, norm=norm)\n    sm.set_array([])  # Empty array, as we'll not use actual data\n\n    grid_size = len(m.data1[0][0])\n\n    plot_one(m.data1[i], 0, cmap, ax)\n    #axes[c].set_title('avg links = {}'.format(self.get_avg_links()[c]))\n    ax.set_xlim(-(1+0.05*grid_size), 2+grid_size*1.05)\n    ax.set_ylim(-(1+0.05*grid_size), 2+grid_size*1.05)\n    \n    #axes[c].axis('square')\n    #plt.show()\n"
In [ ]:
%matplotlib inline
from ipywidgets import interact

# Create an interactive slider from 0 to 99
interact(plot_over, i=(0, 9))
interactive(children=(IntSlider(value=4, description='i', max=9), Output()), _dom_classes=('widget-interact',)…
Out[ ]:
<function __main__.plot_over(i)>
In [ ]:
from matplotlib.animation import FuncAnimation
from IPython.display import display, HTML

def plot_over(i, m, ax):
    normalized = True
    grid_size = len(m.data1[0][0])

    ax.clear()

    for c in range(3):
        num_segments = int(np.max(m.data1[i][c]) + 1) if not normalized else 2
        cmap = create_cmap(c, num_segments)

        norm = Normalize(vmin=0, vmax=num_segments)
        sm = ScalarMappable(cmap=cmap, norm=norm)
        sm.set_array([])

        plot_one(m.data1[i], c, cmap, ax, 0.6)
        ax.set_xlim(-(1 + 0.05 * grid_size), 2 + grid_size * 1.05)
        ax.set_ylim(-(1 + 0.05 * grid_size), 2 + grid_size * 1.05)
        ax.set_title(r'$\beta$ = {}     step {}'.format(1.5, int(i * 10_000))) ###
        ax.axis('off')

fig, ax = plt.subplots(figsize=(12, 12))

# You may need to adjust the range and interval based on your data and desired animation speed
animation = FuncAnimation(fig, lambda i: plot_over(i, m, ax), frames=len(m.data1), interval=250, repeat=False)

# Display the animation in the notebook
HTML(animation.to_jshtml())
Out[ ]:
No description has been provided for this image
No description has been provided for this image
In [ ]:
m.check_state()
Out[ ]:
True
In [ ]:
transformations = [ 
            (1, 1, 1, 1), (-1,-1,-1,-1),      # uniform                        #add later transformation that we can always apply!!!
            (2, 0, 0, 0), (-2, 0, 0, 0),      # single top
            (0, 2, 0, 0), (0, -2, 0, 0),      # single right
            (0, 0, 2, 0), (0, 0, -2, 0),      # single bottom
            (0, 0, 0, 2), (0, 0, 0, -2),      # single left
            (1,-1, 1,-1), (-1, 1,-1, 1),      # swap opposite
            (1, 1,-1,-1), (-1,-1, 1, 1),      # swap adjacent
            (-1, 1,1,-1), ( 1,-1,-1, 1),   
            (1,-1,-1,-1), (-1, 1, 1, 1),      # triple up
            (-1,1,-1,-1), ( 1,-1, 1, 1),      # triple right
            (-1,-1,1,-1), ( 1, 1,-1, 1),      # triple bottom
            (-1,-1,-1,1), ( 1, 1, 1,-1)       # triple left
            ]    

for t in transformations:
    if t[3] < 0:
        print('if {} in transformations: transformations.remove({})'.format(t,t))
if (-1, -1, -1, -1) in transformations: transformations.remove((-1, -1, -1, -1))
if (0, 0, 0, -2) in transformations: transformations.remove((0, 0, 0, -2))
if (1, -1, 1, -1) in transformations: transformations.remove((1, -1, 1, -1))
if (1, 1, -1, -1) in transformations: transformations.remove((1, 1, -1, -1))
if (-1, 1, 1, -1) in transformations: transformations.remove((-1, 1, 1, -1))
if (1, -1, -1, -1) in transformations: transformations.remove((1, -1, -1, -1))
if (-1, 1, -1, -1) in transformations: transformations.remove((-1, 1, -1, -1))
if (-1, -1, 1, -1) in transformations: transformations.remove((-1, -1, 1, -1))
if (1, 1, 1, -1) in transformations: transformations.remove((1, 1, 1, -1))
In [ ]: